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Abstract 

In this latest of our series of Skyrme-HFB mass models, HFB-16, we introduce the new feature 
O ; of requiring that the contact pairing force reproduce at each density the * S pairing gap of neutron 
matter as determined in microscopic calculations with realistic nucleon-nucleon forces. We retain 
the earlier constraints on the Skyrme force of reproducing the energy-density curve of neutron 

>. 

C ! matter, and of having an isoscalar effective mass of 0.8M in symmetric infinite nuclear matter at 
the saturation density; we also keep the recently adopted device of dropping Coulomb exchange. 
Q\ ■ Furthermore, the correction term for the spurious energy of collective motion has a form that 

o. 

is known to favour fission barriers that are in good agreement with experiment. Despite the 



extra constraints on the effective force, we have achieved a better fit to the mass data than any 



^ ■ other mean field model, the rms error on the 2149 measured masses of nuclei with N and Z > 8 
having been reduced to 0.632 MeV; the improvement is particularly striking for the most neutron- 
rich nuclei. Moreover, it turns out that even with no flexibility at all remaining for the pairing 
force, the spectral pairing gaps that we find suggest that level densities in good agreement with 
experiment should be obtained. This new force is thus particularly well-suited for astrophysical 
applications, such as stellar nucleosynthesis and neutron-star crusts. 

PACS numbers: 21.10.Dr,21.30.-x,21.60.Jz,26.60.Gj 
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I. INTRODUCTION 



The r-process of stellar nucleosynthesis is known to depend on the masses and fission 
barriers (among other quantities) of nuclei that are so neutron-rich that there is no hope of 
being able to measure them in the laboratory in the foreseeable future jj. Moreover, the 
composition of the outer crust of neutron stars depends on the masses of nuclei lying as far 
out on the nuclear chart as the neutron drip line [2j. It is thus of the greatest importance 
to be able to make reliable extrapolations of these quantities away from the known region, 
relatively close to the stability line, out towards the neutron drip line. In order to put these 
extrapolations on as rigorous a footing as is feasible at the present time we have constructed 
a series of mass models based on the Hartree-Fock-Bogoliubov (HFB) method with Skyrme 
forces and a contact pairing force {3], 0, O, 3, S ^, ID, U|, the parameters of which are 
fitted to essentially all the available mass data. Clearly, the fit to the mass data must be 
as good as possible; hitherto our best fit was that of model HFB-8 [7|, for which the rms 
deviation with respect to the 2149 measured masses of nuclei with N and Z > 8 given in 



121, is 0.635 



the most recent data compilation, the 2003 Atomic Mass Evaluation (AME) 
MeV. 

However, masses are not the only nuclear quantity of astrophysical interest that can be 
calculated with an effective interaction, and we have been paying an increasing amount of 
attention to such questions as level densities, fission barriers and various properties of the 
inner crust of neutron stars (see Section 5 of Ref. [8j for a brief discussion of the structure 
of neutron stars that is sufficient for our present purposes, and Refs. Q for a more complete 
account). Our ultimate aim is to construct a universal nuclear effective interaction for all 
the various astrophysical applications, and to this end we have been imposing on our mass 
models an increasing number of relevant constraints, even if this entails a slight deterioration 
in the quality of the mass fit. The progress made so far in this respect can be summarized 
as follows. 

i) To make reliable calculations of the inner crust of neutron stars, where neutron-rich 
clusters coexist with a neutron liquid which may be superfluid, it is essential that the value 
of the effective mass implied by the effective interaction be realistic. For its isoscalar com- 
ponent M* the preferred value at the equilibrium density p of symmetric infinite nuclear 
matter (INM) is around 0.8M (see Ref. lOj for a summary of the experimental and theo- 
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retical evidence). While a phenomeno logical value of M*/M much closer to unity had been 
traditionally regarded as essential for good masses in mean-field calculations, we showed {(| 
that it is possible to have fits that are almost as good with M*/M taking the realistic value 
of 0.8. The key to this success lay with exploiting the degree of freedom associated with 
the pairing ; it should nevertheless be stressed that reducing M*/M in this way still has an 
adverse effect on single-particle (s.p.) spectra in the vicinity of the Fermi level p]. 

ii) Beginning with mass model HFB-9 [8] we require that the effective interaction repro- 
duce the energy-density curve of neutron matter, as calculated with realistic 2- and 3-nucleon 
forces [3]. Not only is this condition obviously relevant to all neutron-star applications, but 
it will also improve confidence in finite-nucleus extrapolations out towards the neutron drip 
line. Of all our constraints this was the one having the most adverse effect on the mass fit, 
but it was still possible to hold the deterioration to within tolerable limits. 

iii) By taking a pairing force that was weaker than that of model HFB-9 and of all earlier 
models we found that it was possible, for the first time with our models, to obtain level 
densities in reasonable agreement with experiment {9], [jjj] , while maintaining an acceptable 
mass fit. 

iv) By making a phenomenological adjustment of the collective correction it is possible 
to fit both masses and fission barriers within the same model (HFB-14) {lo| . 

v) Our latest model (HFB-15) |TJj incorporates the frequently recommended procedure 
of dropping Coulomb exchange, thereby simulating neglected effects such as Coulomb cor- 
relations, charge-symmetry breaking of the nuclear forces, vacuum polarization, etc. A 
significant improvement in the mass-data fit was achieved. 

As an alternative to non-relativistic Skyrme-based mass models, mass models based on 
the relativistic mean-field (RMF) method have also been proposed Q] ■ Now the isospin 
dependence of the spin-orbit field is formally quite different in the two approaches, and 
it might be thought that RMF models would provide more reliable extrapolations to the 
neutron drip line, since one must inevitably prefer a model that respects Lorentz covariance 
to one that does not, other things being equal. However, it has been shown that if the same 
mass data are fitted equally well in the two approaches then essentially the same predictions 
are made out at the neutron drip line [nj]. Thus we believe that the Skyrme approach is just 
as reliable as the RMF method; in any case, the latter method has not yet been developed 
to the point where its mass fits are as good as those given by Skyrme-based methods. 
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In the present paper we impose the additional constraint force of requiring that the 
effective interaction reproduce as a function of density the microscopic 1 S'o pairing gap of 
neutron matter, while retaining all the constraints of the earlier papers. Without this new 
constraint it would be impossible to make reliable investigations of a possible superffuid 
phase in the inner crust of neutron stars. For this purpose, we develop a scheme first 



proposed by Duguet [18| for constructing an effective density-dependent contact pairing 
interaction which reproduces any given microscopic gap exactly. Instead of parametrizing 
the density dependence by a simple functional form, the strength of the contact interaction 
is determined at each density by solving directly the HFB equations in neutron matter, 
and requiring that the resulting gap be equal to the given microscopic gap at the Fermi 
momentum corresponding to that density. 

In Section [IT] we describe the way in which we implement the HFB method, with particular 
emphasis on points that we overlooked in our earlier description in Ref. [2]]. Section UTT1 
explains our procedure for determining the strength of the pairing force, while the results of 
the new mass fit are presented in Section IIV1 The two appendices summarize some relevant 
features of the HFB method. 



II. THE HFB METHOD 

The Skyrme force that we use in the present calculations has, as with all our previous 
HFB mass models, the conventional form 

v^farj) = t„(l + x P a )5( rij ) + + x 1 P a )^ [pi 8{r i5 ) + <$(»"<,•) 
-K 2 (l + x 2 P (T )^2Pij-S(r i j)p ij + h 3 (l + x 3 P a )p{ry 5{r i5 ) 
+^Wo((Ti + (Tj) -Pij x 5{r ij )p ij , (2.1) 

where rij = Vi —rj, r = (r^ +rj)/2, = — ifr(Vi — Vj)/2 is the relative momentum, and P a 
is the two-body spin-exchange operator. The contact pairing force that we take here acts, as 
before, only between nucleons of the same charge state q (q = n or p for neutron or proton, 
respectively) 

v^(r i ,r j ) = v^[p n (r),p p (r)]5(r lJ ) , (2.2) 
where v nQ [p n ,p p ] is a functional of the nucleon densities that will be specified in Section [TTT] 
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The only paper in which we attempt to describe the way in which we implement the 
HFB method is Ref. [3] (Section 2). Since this account was somewhat incomplete, and in a 
sense misleading, we present here a more detailed description of what we actually did. Our 
original presentation [2I was developed within the framework of the standard formulation 



of Mang [19| and of Ring and Schuck 20], which is certainly well adapted to the use of a 
discrete basis such as an oscillator basis, as in the case of our own calculations. Nevertheless, 
some clarification is necessary when this formalism is applied to effective forces of the form 
given in Eqs. (J2ITJ and (12T2D . 



The standard formulation 
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201 ] assumes two-body forces, and starts with the expression 



H = l v c \ c 3 + 4 ^2 v HM c \ c ) c i c k (2-3) 

ij ij,kl 

for the Hamiltonian (see, for example, Eq. (5.25) of Ref. [3]). Here we are working in a 
fixed basis of discrete s.p. states labelled by i,j, etc., e.g., an oscillator basis, with c\(cj) 
denoting creation (destruction) operators for real nucleons in such states (the charge type is 
implicit in the label). Also we have introduced the antisymmetrized matrix element of the 
two-body force 

Vij,kl = Vij,kl ~ v ij,lk , (2.4) 

and the matrix element t%j of the kinetic-energy operator — fr 2 V 2 /2M g , denoting the nucleon 
mass by M q . 

The Bogoliubov transformation 



I 

then defines creation (destruction) operators Pl(Pk) of quasiparticles as linear combinations 
of the creation and destruction operators of real nucleons. The essence of the HFB method 
is that it takes the vacuum state of these quasiparticles as the trial function for a variational 
approximation to the energy of the nucleus. Denoting this (normalised) vacuum state by 
1^ >, and noting that it will be a function of the transformation coefficients Uik , V%j , etc., 



the HFB ground-state energy becomes, as for example in Eq. (E.20) of Ref. 20j 



E mu =< *\ fl\* > ■ Tr ( tp + il> - iAK* j . (2.G) 



in which Tr denotes the trace, p and k are the so-called normal and abnormal density 
matrices, respectively, 

Pij =< tflctcil* >= V u V H = P% (2-7a) 



and 



while 



Kij =< *|c,-Ci|* >= ^ VJE/j, = /Cjj , (2.7b) 



rjw = 22 v ki M Pii (2.8a) 



and 



A w = - ^ /% . (2.8b) 

ij 

We now notice the first departure from this standard formulation that we have to make on 
account of our choice of effective forces of the form given in Eqs. (12.11) and (12. 2p . Although 
Eq. (12.61) remains valid, along with Eqs. (12 .4p . (I2.5p . (j2.7aj) and (12.7bj) . because the force in 
the particle-particle channel is different from the one used in the particle-hole channel we 
must replace Eqs. (I2.8al) and ( I2.8bl) . respectively, by 

r« = £e2^' + E e 2#ft-'' (2.9a) 



where v^f) are the antisymmetrized matrix elements of the Coulomb interaction, and 



^EC^ > ( 2 - 9b ) 



ki 2 

in which the choice of force in the respective channels is made explicit. This distinction was 
made clear in our first HFB paper [3j (see Eq. (7) of that paper, where, however, there was 
an error in the subscripts of the expression for Vk P ), but it was incorrect to imply that we 
could still take a Hamiltonian of the form (12.31) as the starting point: this would be possible 
only if no matrix element Vki,ij appeared in both Fj-i and A^, which is not the case. Rather, 
one should really proceed via a density-functional approach, as discussed below. 
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The second departure from the standard formulation that we are obliged to make because 
of our choice of effective forces was not made clear at all in Ref. The HFB equations 
resulting from the minimization of the total ground-state energy (12. 6p are 

y ( ^ ~ X6ij Atj ) ( Ujk ) = E J Uik ) (2 io) 

i V ~ A h -Kl + Ma) \v* ) \v ik ) 

where E are the quasi-particle energies, and A is the usual chemical potential, arising as a 
Lagrange multiplier. Also 

K i = 9 -^ = h' h (2.11a) 

denoting matrix elements of the self-consistent s.p. Hamiltonian, and 

A y = ^^ = -A„, (2.11b) 

denoting matrix elements of the pairing potential. Now Eq. (12.101) differs from Eq. (6) of 
Ref. js} in that the quantity hij defined in the first member of Eq. (7) of that paper, 



tij ~t~ 7 (2.1 2) 



has been replaced here by 



hl^Kj + h]™, (2.13) 



where we have introduced matrix elements of the rearrangement s.p. field 

h rear _ 1 ( UU kl,pm 1 UU kl,pm * (r> -, ^ 

=2^[^~ PmlPpk ~2^~ KpmKlk ' (2 - 14) 

klpm \ J J / 

which arises entirely from the derivative with respect to density of the density-dependent 
components of the force (the Coulomb interaction is independent of the density and therefore 
does not lead to rearrangement terms). Although these rearrangement terms were omitted 
from Eq. (6) of Ref. [3j , they were included in all our actual calculations. (The rear rang ement 



terms were also dropped from Eq. (7.42) of Ref. [20], and from Eq. (14) of Ref. 58] . Note 
also that Eq. (7.42) of Ref. 3 

absorbs A into h.) Let us remark that the forces considered 
in this work depend only on the normal nucleon densities but not on the abnormal densities, 
whence the expressions (I2.11bl) and ( I2.9bl) are equivalent. The foregoing formalism could 
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certainly be implemented if one calculated the two-body matrix elements v^ y kl and w? a ^ for 
forces (12. II) and (12.2h . respectively, in the oscillator basis. However, although this presents 
no problem of principle, it is not what one does, despite suggestions to the contrary in Ref. 
{sj]. Rather, one follows a simpler approach that does not require the calculation of these 
matrix elements, even if one still solves the HFB equations (12.101) in the oscillator basis, 
as we do here. Instead of starting with the Hamiltonian (I2.3p . one begins by noting that 
for zero-range interactions, such as forces (12 .ip and (12. 2p considered in this work, the HFB 
energy can be written as the integral of a purely local energy-density functional 



#hfb = J £ H FB(r)d d r , (2.15) 
where, assuming time reversibility, as always in this paper, 

£hfb(>) = £sk y Pn{r), Vp n {r),r n (r), J„(r),p p (r), V p p {r) , r p (r) , J p (r) 



+ £, 



Coul 



p p (r) +£ P air Pn(r),p n (r),p p (r),p p (r) , (2.16) 



in which the usual normal density p g , kinetic-energy density T q , spin-current density J q and 
pairing density p q appear (see Appendix [A] for precise definitions). 

Expressions for the first term in Eq. (I2.16p . the energy density associated with the Skyrme 



force (12.11) and implicitly including the kinetic energy, can be found, for example, in Ref. 2l|. 
The second term in Eq. (12.161) represents the Coulomb-energy density, the exchange part of 
which is treated in our code by the Kohn-Sham approximation [22| (usually referred to as 



the Slater approximation in the nuclear-physics literature, but the original Slater version [23] 
is somewhat different), assuring thereby that the energy- density functional remains local. 
Actually, in the model described in the present paper, HFB-16, we drop the Coulomb- 
exchange term, as discussed in Section 0], and thus write 

£coui = ^e Pch V Coul , (2.17) 

in which ep c h is the charge density associated with protons (this differs from ep p because we 
are taking account of the finite size of the proton), and V Conl is the electrostatic potential, 
given by 

j 



V c ™\r) =e / (Ft>£*£-L. (2.18) 



r — r 
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The last term in Eq. (12.161) is the pairing-energy density, associated with the pairing 
force (12. 2p . and is given by 

= \ Yl v* q [ P n(r), P P (r)]p q (r) 2 . (2.19) 

q=n,p 

By taking a density functional of the general form given by Eqs. (12.161) . rather than 
the Hamiltonian (12.31) . as the starting point, the possibility of having different equivalent 
forces in the particle-hole and particle-particle channels arises quite naturally. Indeed, in 
this formulation the forces never have to appear explicitly, and one even has the freedom 
to choose the form of the functional in such a way that no corresponding forces exist at all. 
However, in this paper our choice of functional does correspond rigorously to forces of the 
form ( 12.11) and ( 12.21) . as we have already indicated. 

When the HFB energy is expressed in terms of an energy- density functional, as in 
Eqs. (12.151) and (I2.16p . then minimizing it with respect to the normal and pairing density 
matrices, for a fixed average number of neutrons and protons, leads to the HFB equations 
in coordinate space, as shown by Dobaczewski and co-workers 24], |25( (the link with the 
discrete-basis formulation (12.101) is discussed at the end of Appendix |A|) . The advantage of 



solving the HFB equations in coordinate space, as in Refs. 24j, |25j, is that it facilitates an 
accurate determination of the asymptotic wavef unctions, but it has little impact on the de- 
termination of binding energies, which are our main concern. Thus for convenience we solve 
the HFB equations in the discrete-basis form (I2.10p . An essential step in this solution is the 
computation of the oscillator matrix elements (12. llaj) and (12.11bj) . but we considerably sim- 
plify the numerical task by generalizing to HFB the procedure outlined by Vautherin for the 
first HF-BCS calculation of deformed nuclei 26J]. This procedure (see Appendix [A] for more 
details) works essentially in an oscillator basis but makes a detour into the coordinate-space 
representation by expressing these matrix elements as matrix elements of the self-consistent 
s.p Hamiltonian h ' [r) a i a and the self-consistent pairing field A 9 (r), respectively, both of 
which appear in the coordinate-space form of the HFB equations ( 1A.18I) . (explicit expres- 
sions for both these fields are given in Appendix |A|). In this way the computation of the 
far more numerous matrix elements H2.4I) of the two-body force is no longer required. From 
this point on, our calculations proceed on the general lines described in Ref. jsj. 
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III. DETERMINATION OF DENSITY DEPENDENCE OF PAIRING FORCE 



Little is known about the density dependence of the effective pairing interaction. The sim- 
plest ansatz, commonly referred as "volume pairing", is simply to assume that v 7Tq [p n ,p p ] = 
V 1Tq ) a constant, independent of the density, in which case the corresponding rearrange- 
ment term vanishes. This prescription has been adopted for all our previous forces except 
for the parameter sets BSk3 5], BSk5 (| and BSk7 [6], for which we considered a density 
dependence of the usual form 



V 7rQ [p n ,p p ] = V nq <l-rj 



(3.1) 



where p = p n + p p , while V nq , rj and a are adjustable parameters. Testing different sets of 
parameters, it was concluded that this particular form (13. 1 p does not improve the global fit to 
the experimental mass data, and we thus abandoned it in our subsequent parametrizations. 
However, if one is to reproduce not only masses but also the microscopic pairing gaps given 
by bare nucleon-nucleon potentials, i.e., interactions fitted directly to the 2- and 3-nucleon 
data, then the effective pairing force must be density dependent, as was first shown by 
Garrido et al. {27] with a contact pairing force with the functional form (13.11) . 

In the present paper, we wish to constrain our effective force to the pairing properties of 
pure neutron matter. Several microscopic many-body calculations of the 1 Sq pairing gaps 
have been carried out using various approaches and approximation schemes. Calculations at 
the simplest BCS level with the free s.p. spectrum for different nucleon-nucleon potentials 



yield very similar results 



28 



J, [30|. In fact, at this level the pairing gap is essentially 
determined by the experimental 1 S phase shifts 3l| • Calculations going beyond the BCS 
approximation generally show that medium effects lower the maximum of the gap function 
A n (p n ) but they disagree among themselves as to the precise density dependence. (Note 
that we use the notation A n (p n ) for the pairing gap evaluated at the Fermi momentum 
f^Fn — (3vr 2 p„) 1 / 3 .) Given this lack of agreement, we take here for the microscopic pairing 



gap the one shown in Fig. 7 of Ref. 



given by the realistic Argonne v 



"orce in the 



32 



33 



M 



35] 



BCS approximation. In any case, recent quantum Monte-Carlo calculations 
suggest that the "true" gap might be actually close to the BCS one. 

We show this microscopic BCS g function of neutron-matter density in Fig. [TJ for 
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FIG. 1: So pairing gap A n of neutron matter as a function of the density p n calculated at the 
Fermi momentum kp n = (Svr 2 /^) 1 / 3 for the realistic Argonne V14 interaction and various 

effective forces. Note that the effective force BSkl6 derived in this paper yields a pairing gap 
indistinguishable from the curve labelled "realistic V14 (BSkl6)". 



convenience, we have used the essentially exact analytical representation 36] 



A„ W = - fe.) A.^^J^M. , (3. 2) 

where 9 is the Heaviside unit-step function, and the associated parameters are given in 
Table [H For comparison, we also display in the same figure the gaps corresponding to 
various effective interactions. The large differences between the neutron matter gaps for our 
forces BSk7 and BSk8 is particularly striking, given that they both yield good mass fits and 
hence comparable pairing gaps in finite nuclei. The gap obtained with the force BSk7 is 
closer to the microscopic one, which might be related to the fact that its pairing component 
is density dependent, but the discrepancy with respect to the microscopic gap is still quite 



unacceptable. Much better agreement is found for the Skyrme force SLy4 37|], using the 



pairing adopted in Ref . [38j . However, the rms deviation of the mass fit given by this model 



is very large (5.1 MeV for the subset of even-even nuclei 38|, |39(). In the present paper 
we show that despite the failure of the BSk7 and BSk8 models to yield acceptable pairing 
gaps in neutron matter, it is still possible to construct a modified BSk model that fits the 
microscopic pairing gap of neutron matter exactly, while maintaining the quality of the mass 
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fits obtained with our earlier models. We label the new model BSkl6, as indicated in Fig. [TJ 
Now if we adopted the parametrization (13 .ip for the density dependence of the contact 
pairing force in our new mass model a serious problem would arise. In any calculation of 
this sort the values of the pairing parameters depend on the effective interaction adopted in 
the particle-hole channel, essentially through the effective mass (see Eqns. (13. 3p and (13.41) 
below). In the case of the work of Garrido et al. j^, a specific force (Gogny) was chosen at 
the outset as the interaction for the particle-hole channel, whence the parameters V^q, rj and 
a of the pairing force (13.11) could be easily fitted to a microscopic pairing gap. However in 
the present work, the parameters of the Skyrme force acting in the particle-hole channel are 
not known a priori but rather are determined a posteriori by fitting to experimental nuclear 
mass data, as discussed in Section HVl Fitting the parameters of the effective interaction in 
both the particle-hole and particle-particle channels, while simultaneously reproducing the 
microscopic pairing gap A q (p q ), would be an extremely onerous numerical task when using 
a phenomenological functional, such as Eq. ( 13. ip . Moreover, there is no guarantee that the 
parametrization (13. ip of the density dependence will be optimal in general, as indicated in 
Ref . 4(3] , even though it works well in the case considered by Garrido et al. 27|. Both of the 
above problems are avoided in the fitting procedure that we adopted in the present paper, 
and now describe. 

Following Duguet [3] and more recent calculations 41], we assume that the pairing force 



for nucleons q depends only on p q , and not on the total density p. Besides, we suppose that 
the effective interaction v 7Tq [p q (r)} at the point r is the same as that in INM of pure nucleon 
species q at the density p q (r). But instead of postulating a simple functional form for the 
density dependence, we have determined the strength of the effective pairing force v 7rq [p q ] 
at each nucleon density p q by solving the HFB equations in uniform matter and requiring 
that the resulting gap reproduce exactly the given microscopic pairing gap A q (p q ) at that 
density. The only free parameter is the energy cutoff e\. Note that in the present work, 
we have set an upper limit on the s.p. states that are retained in the spectrum without 
using the Bulgac-Yu [421 ] regularization procedure, which we adopted in models HFB-12 to 
HFB-15 (2,0 

In uniform matter, the HFB equations reduce to the BCS equations (1B.8h and (1B.9[) . 
Solving these equations for v 7Tq [p q ] in pure matter of nucleon species q and adopting the 
cutoff prescription (iii) discussed in Appendix [HI yields (after a change of variable £ = e — U q 
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in the integrals) 

v **\p q ] = -8tt 2 ( ^ ' ( r +eA d£-== ^ \ , (3.3) 

where M*(p q ) is the effective nucleon mass in matter of pure nucleon species q at density 
p q , as given by Eq. (IA.14D . and p q = \ q — U q can be obtained by solving 



Note that determining the strength f 719 of the effective pairing force by Eq. (13.31) ensures its 
automatic renormalization for any changes of the energy cutoff £\. The same value of E\ is 
used in calculations of finite nuclei (see Eq. (14. 3 j) below). Let us remark that if, instead of 
imposing a cutoff above the Fermi level, we would have chosen a fixed cutoff (as was done 
in model HFB-1 [3!]), Eq. (13. 3p would have been actually an integro-differential equation for 
■u 71 " 9 !/^] as discussed in Appendix El In such case, it would therefore have been much more 
difficult to compute the strength of the effective pairing force. 

In contrast to the use of parametrized pairing forces (13.11) . in the present scheme we 
reproduce exactly the density dependence of the microscopic pairing gap. At the same time, 
the fitting procedure of the Skyrme parameters is greatly simplified: at each iteration of 
the mass fit, the pairing function v nq [p q ] is determined unambiguously through the coupled 
equations (13.31) and ( 13.41) by the given microscopic gap A n (p n ) and the running values of the 
Skyrme parameters ti, t 2 , X\ and x 2 . 

In our actual numerical calculations, instead of solving both Eqs. (13. 3p and (13.41) at each 
density, we approximate the reduced chemical potential p q by the Fermi energy 

(g) _ ^Fg ( 5) 

- 2M* ' [6 - b) 

where kp q is the Fermi wave number 

k Fg = (^ 2 p q ) 1/3 ■ (3.6) 

We have found that this approximation holds well, provided A q ^ £p . Then we do not 
have to solve Eq. ( 13.4D at all, and the evaluation of v nq [p q ] through Eq. ( I3.3P reduces to 
a single, simple integration. The effective pairing strength v 7Tq [p q ] determined for the final 
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FIG. 2: Density dependence of the strength of the effective pairing force v 7TQ [p q ], which has been 
fitted to the 1 Sq pairing gap A n shown in Fig. [TJ The effective mass at each density is calculated 
with the final BSkl6 Skyrme parameters (see Section ITV|) . For comparison, we show the effective 



pairing force used in the SLy4 mass model 



38, 



39]. 



Skyrme force BSkl6 (see Section HVl) through Eq. (13. 3p . is shown in Fig. [2] it will be seen 
that the density dependence is qualitatively similar to that given by Eq. (13.11) . 

The foregoing neutron-matter prescription suffices in principle to fix the neutron pairing 
(see, however, Section [IV] for the distinction that we make between even-iV and odd-iV 
nuclei). Because of charge-symmetry breaking, the proton pairing force should in principle 
be determined in an analogous way on the basis of "proton-matter" calculations, which are 
performed in exactly the same way as are neutron-matter calculations, using the appropriate 
proton-proton interaction, with, of course, the Coulomb force suppressed, since "proton- 
matter" calculations would otherwise diverge. However, this prescription will not suffice, 
since in finite nuclei the Coulomb force does act, and may be expected to modify the proton 
pairing, in a highly complicated way. For simplicity, we take v np [pp\ to be given by v 7Tn [p p ], 
multiplied by a constant, density-independent factor which is the same for all nuclei and is 
taken as a fitting parameter (see Section [TV] for more details). 
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FIG. 3: Deformational variation of collective correction of Eq. (|4.2[) . 
IV. THE MASS FITS 



We have described our treatment of the Skyrme and pairing forces of our new model, 
HFB-16, in Sections HT1 and HTT| respectively, but, before presenting the parameter sets for 
these forces that emerge from the fits to the mass data, we discuss several further points, as 
follows. 

i) Wigner correction. To the HFB energy calculated for the Skyrme force (12. ip and the 
pairing force ( 12. 2ft has to be added the Wigner correction 



E w = V w exp 




(4.1] 



which contributes significantly only for light nuclei [4|. Our treatment of this correction 



is purel y p henomenological, with the first term believed to be representing a T = n — p 
pairing 



43j, |44], |45j , while the second is characteristic of Wigner's supermultiplet theory |46j, 



based on SU(4) spin-isospin symmetry, although it can also be interpreted as corresponding 
to T = 1 n — p pairing 43] . 

ii) Collective correction. As in all our previous models we subtract from the calculated 
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HFB energy an estimate for the spurious collective energy. The form we adopt here for this 
correction is 

E coll = EZt nk {b tanh(c|/3 2 |) + d\fa\ exp{-l(\fa\ - & ) 2 }} , (4.2) 

in which E^ nk denotes the cranking-model value of the rotational correction and fa the 
quadrupole deformation, while all other parameters are free fitting parameters. This cor- 
rection is intended to take account of both rotational and vibrational spurious energy, but 
while this expression for E co u is seen to vanish for spherical nuclei, we know that a vibra- 
tional correction is still required for such nuclei. Thus we must suppose that the vibrational 
correction for spherical nuclei is absorbed into the fitted force parameters, so that it is only 
the deformational variation of the vibrational correction that is represented by Eq. (14. 2p . 
along with the complete rotational correction. Our final mass fit gives for the coefficient b 
appearing in the first term of Eq. (14 .21) the value of 0.8 (see Table ITTTT) . the same value as we 
found for model HFB-14 [h]]. Thus the argument that was made in Ref. 10] points once 



again to the first term of Eq. (14. 2p being identified with the rotational correction. 



Actually, the second term of the correction (14. 2p differs from that of model HFB-14 [lOj in 
that it carries an additional factor \fa\, the role of which is to ensure that this term, like the 
first term, vanishes as sphericity is approached. Some instabilities that were encountered 



for spherical or weakly deformed nuclei with HFB-14 [lOj are thereby avoided, and this 
modification is in part reponsible for the improved mass fit of this new model (see below in 
this Section). The corrections for the two models are compared graphically in Fig. [3l where 
it will be seen that except in the spherical limit the differences between the two corrections 
are small, especially for large deformations. This point is of crucial importance for fission 
barriers, which are highly sensitive to the collective correction at large deformations, and we 
do not expect that the barriers to be calculated (in a future paper) with the present model 
will be appreciably different from those found with model HFB-14. In any case, it will be 
possible to perform any fine tuning of the barriers that may be required by adjustment of the 
collective correction alone, without any modification of the Skyrme or pairing parameters, 
our experience in Ref. |ljj] showing that the perturbation of the mass fit will be minimal, 



provided the collective contribution has the form ( 14.2 



iii) As with models HFB-14 



10| and HFB-15 [11), we impose on the parameters of the 



Skyrme force the condition M*/M = 0.8 for the isoscalar effective mass in symmetric INM 

16 



at the density p = po (see Section [IJ . 

iv) By varying the exponent 7 the incompressibility coefficient K v can be varied, and we 



find that excellent mass fits can be obtained over the range 230 < K v < 270 MeV. £ 
measurements of the breathing mode restrict this range to 230 < K v < 250 MeV 
so we limit the value of 7 accordingly. 



owcvcr, 



47 



48|, 



v) The measured rms charge radius R c of 208 Pb, 5.501 ±0.001 fm 49], was required to 



be well reproduced; relaxing this condition leads to very little improvement in the mass fit. 



vi) As with all our models since HFB-9 



we constrain our Skyrme force to fit neutron 



matter (see Section [II); this turns out to be equivalent to setting the INM parameter J to 
be 30 MeV. 

vii) We drop the Coulomb-exchange term from the HFB energy-density functional (see 
Section [I]). 

viii) As with all our models, we treat the case of odd N and/or odd Z by the "level- 
filling" approximation. That is, the odd nucleon is placed with equal probability in each of 
the lowest-energy available degenerate states generated by the HF calculation with the next 
lowest even number of nucleons, and then applying blocking. 

ix) Our new prescription for the pairing strength, described in Section [Till is strictly valid 
only for neutrons. Because of Coulomb effects, and a possible charge-symmetry breaking of 
nuclear forces, we must allow for the proton pairing strength to be different. Likewise, we 
shall allow the pairing to be different for odd-A and odd-odd nuclei to compensate for our 
failure for such nuclei to project out states that respect time-reversal invariance. We take 
account of these extra degrees of freedom by multiplying the value of v wq [p q ], as determined 
through Eq. fl3.3|) . by renormalizing factors , where fp,fp and /~ are free, density- 
independent parameters to be included in the mass fit; in keeping with the spirit of this 
paper we set /+ = 1. 

x) In calculations of finite nuclei, we use the smooth pairing cutoff factor 

-1/2 



ft 



1 + exp((e l - A-e A )/r) 



(4.3) 



with r = 0.25 MeV, in order to prevent an unphysical selection of nearly degenerate levels 
whenever the cutoff energy lies between the levels. 

xi) The spurious centre-of-mass energy is removed following the essentially exact proce- 
dure described in Ref. Q]. 
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xii) A correction for the finite size of the proton is made to both the charge radius and the 
energy, as in all our previous HFB models. We assume a Gauss distribution of charge over 
the proton, with an rms radius of 0.895 fm 50j; the folding of this charge distribution over 
the HFB distribution of point protons is performed using Eqs. (4.2) and (4.3) of Ref. 

xiii) As in all our papers except Ref. [7| we make no attempt to project out states of 
good particle number. 

xiv) The oscillator basis in terms of which the HFB wavefunctions are expanded contains 
20 major shells. 

The new force parameters, labelled BSkl6, are fitted to the same set of 2149 measured 

n n 

masses [12[ as are all our models since HFB-9 [8j] ; the resulting values are shown in Table [TT1 
while Table IIHI shows the parameters of the collective correction of Eq. (14. 2p . These two 
tables define the HFB- 16 model, with which we have constructed a complete mass table 
running from one drip line to the other over the range Z and N > 8 and Z < 110. 

The first ten lines of Table HT1 show the Skyrme parameters; their only noteworthy feature 
is the small value of and the large value of X2, which implies that the X P and 3 P interactions 
have strengths of comparable magnitude but opposite sign, the latter being attractive. The 
next three lines show the factors by which the pairing strength determined from neutron 
matter must be renormalized for protons and for odd nuclei. With /+ = 1 we see that to 
within the limits of numerical accuracy f^/fn = fp/fp an d the proton and neutron pairing 
strengths are effectively equal. This means that these three degrees of freedom could have 
been reduced to a single one, the ratio of the odd-number pairing strength to the even- 
number pairing strength (as in all our previous models, pairing is always a little stronger 
for an odd number of nucleons). Line 14 shows the pairing cutoff parameter E\; the value of 
16 MeV was chosen on the basis of past experience {{J. The last four lines give the Wigner 
parameters of Eq. (14.11) . 

The rms and mean (data - theory) values of the deviations between the measured masses 
and the HFB-16 predictions are given in the first and second lines, respectively, of Table [TVT. 
where we also compare with HFB-15 [11], HFB-14 10J and our previous "best-fit" model 
HFB-8 We see that with this new model we have achieved our best fit ever (er = 0.632 
MeV), in addition to having a greater conformity to physical reality than with any of our 
other models. It should be recalled that the first physical constraint we imposed on our mass 
fits, conformity to the neutron- matter energy- density curve, led, with model HFB-9 js], to 
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a serious loss in the quality of the mass fit. However, we have now entirely recovered the 
precision of model HFB-8; this could be a result of the device of dropping Coulomb exchange 
(introduced in model HFB-15 [11]), of modifying the collective correction discussed above, 
or of requiring a better agreement with the BCS neutron-matter pairing gap. Concerning 
this last point, we should emphasize that starting from a neutron- matter gap that includes 
medium effects beyond the BCS approximation, we were not able to obtain a good mass fit. 
Fitting masses to such microscopic gaps seems to be very difficult in the present framework. 

The superiority of the new model becomes still more striking on looking at the next six 
lines of Table IIVI Lines 3 and 4 show the rms and mean deviations for the astrophysically 
crucial subset of the mass data consisting of the 185 neutron-rich nuclei having a neutron- 
separation energy S n < 5.0 MeV. The following four lines give the rms and mean deviations 
for the S n and the beta-decay energies Qp, using the full data set of 2149 measured masses: 
since these are differential quantities they are astrophysically more relevant than the absolute 
masses M. It is noteworthy that in all these important categories of mass-related quantities 
model HFB-16 out-performs all our other HFB models. 

Turning to charge radii, lines 9 and 10 show the rms and mean deviations between the 



measured values 



49] and the model predictions. Only here do we see any deterioration 



with respect to model HFB-8; the mean deviation e(R c ) suggests that from model HFB-14 
onwards we should have been taking a slightly larger density po, but we do not know whether 
this would have an adverse effect on the mass fit. 

As stated in Section [J, one of our concerns has been to have a pairing force that not only 
permits a good mass fit but also is consistent with level densities that are in reasonable 
agreement with experiment. However, in the present model the pairing force is fixed almost 
entirely by the a priori calculated pairing properties of neutron matter, the flexibility offered 
by the —Nation fetors f* being neghgib.e in this respeet. Now ptev.ons experience 
with model HFB-13 [9|, H4| shows that reasonable level densities are found when the spectral 
pairing gap (uvA) of the model (defined as in Eq. (5) of Ref. js]) lies close to the fifth-order 
experimental even-odd mass differences A^ 5 ^. We compare in Fig. H] the HFB-16 spectral 
pairing gaps for the Sn and Pb isotope chains with those of model HFB-13 [9]. It will be 
seen that the new model HFB-16 agrees closely with the earlier one in the case of the Pb 
isotopes, and is not substantially stronger in the case of the Sn isotopes (similar conclusions 
hold for the third-order even-odd differences A^ 3 '). We leave for a later paper the actual 
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FIG. 4: Comparison of the experimental even-odd differences A^ 5 ** with the HFB-16 theoretical 
neutron spectral pairing gaps (uvA) for the Sn and Pb isotopic chains. 



level-density calculations with HFB-16, but Fig. @] is encouraging. 

Magic neutron- shell gaps. While the global fit to the mass data given by model HFB-16 
is seen to be excellent, a few individual predictions could still be quite anomalous without 
having a serious impact on the global fit. Of particular importance in this respect are the 
masses involved in the definition of the neutron-shell gaps, 



A n (JV , Z) = S 2n (N , Z) - S 2n (N + 2, Z) 



(4.4) 



In Figs. [5]- [8] we show these gaps as a function of Z for the magic numbers iVo = 50, 82, 126 
and 184. The agreement with experiment is excellent for iVo = 50 and 82, with strong gap 
quenching predicted as the neutron drip line is approached. On the other hand, for iVo = 126 
the agreement with experiment is quite bad, as is always the case with our HFB models. 
No quenching is predicted for either iVo = 126 or 184. 

Macroscopic properties. Table IVl shows the parameters of infinite and semi-infinite nuclear 



matter (INM and SINM) that we have ca 
of these parameters see, for example, Ref. 



culated for the force BSkl6 (for the definition 
j9|). Our SINM calculations, which do not take 



20 




FIG. 5: Nq = 50 shell gap as function of Z for mass model HFB-16. 




FIG. 6: Nq = 82 shell gap as function of Z for mass model HFB-16. 
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FIG. 7: iVo = 126 shell gap as function of Z for mass model HFB-16. 
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FIG. 8: iVo = 184 shell gap as function of Z for mass model HFB-16. 
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pairing into account, are performed using the HF code of M. Farine, as described in Appendix 



A of Ref. 



52j | . It should be noted that the values of M* and J were imposed, as described 



at the beginning of this Section. In this table we show also the corresponding parameters 
for our last two forces; it will be seen that there has been very little change, although it 
is interesting to note that with the omission of Coulomb exchange (BSkl5 and BSkl6) the 
binding energy per nucleon of symmetric INM, a v , has increased. All our values of these 
parameters are compatible with all the available experimental data, as discussed in Ref. 9J. 
To this earlier discussion we would like to append the remark that with all our forces we find 
an isovector effective mass M* that is smaller than the isoscalar effective mass M* at the 
density po- This result implies that the neutron effective mass M* is larger than the proton 
effective mass M* in neutron-rich matter. Such an isovector splitting of the effective mass 
is consistent with measurements of isovector giant resonances 531 ]. and has been confirmed 



in several many-body calculations with realistic forces 5J, [55]. 



The energy-density curve of neutron matter for force BSkl6 is indistinguishable from the 
realistic curve of Ref. [k^] up to the supernuclear density of 0.3 neutron. fm -3 , as is the case 
with all our other forces that have been fitted to J = 30 MeV, i.e, BSk9 and all later forces 
(see Fig. 13 of Ref. |9]). It is to be noted that unlike Ref. [53j we have not had to resort 
to a second £3 term in the Skyrme force in order to simultaneously fit neutron matter and 
obtain the correct sign for the isovector splitting of the effective mass. 

Fig. [9] shows the potential energy per particle in each of the four two-body spin-isospin 
(S, T) channels as a function of density for symmetric INM; we give results for both BSkl6 
and Brueckner-Hartree-Fock (BHF) calculations with realistic two- and three-nucleon forces. 
These latter results are taken from Fig. 6 of Ref. 53J and quoted as Ref. 62 of that paper. 
There is reasonable agreement between BSkl6 and the realistic calculations in all states 
except the (1,1) state, where BSkl6 is strongly attractive, while BHF is very weak. It seems 
to be very difficult, in fact, to fit both odd states with a conventional Skyrme force. For 
example, SLy4 [37| handles the (1,1) state rather well, having imposed %2 = -1, but as a 
result works much worse than BSkl6 for the (0,0) state. In any case, the excessive attraction 
that BSkl6 gives in the (1,1) state leads to the onset of a ferromagnetic instability in neutron 
matter at relatively low density, as seen from the value of Pf rm g/Po in Table PVT. 

Stability of extrapolation. Although all our models agree fairly closely in the known region 
of the nuclear chart, simply by virtue of their fit to the data, there is no a priori guarantee 
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FIG. 9: Potential energy per particle in each (S,T) channel as a function of density for symmetric 
INM. The full symbols correspond to BHF calculations and the open symbols (connected with solid 
lines) to the BSkl6 force. 

that they will give comparable extrapolations to the astrophysically interesting region close 
to the neutron drip line. Accordingly, in Table I VI I we compare the predictions made by 
HFB-16 with our earlier models HFB-15, HFB-14 and HFB-8 for all those nuclei with 
26 < Z < 110 for which S n < 4.0 MeV. This table shows that for the more astrophysically 
relevant quantities S n and Qp the differences between the various model predictions are 
much smaller than for the absolute masses, and in fact are comparable to the deviations 
between each model and the data. This stability of the predictions against slight changes 
in the model is consistent with an overall reliability of the HFB approach, although insofar 
as differences between the model predictions are significant, we would prefer HFB-16, given 
both the excellence of its fit to the data and the high degree of reality built into the model. In 
this latter respect, the conformity of both the Skyrme and pairing components of the effective 
force to realistic neutron-matter calculations should particularly enhance our confidence in 
predictions for highly neutron-rich nuclei. 
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V. CONCLUSIONS 



The new feature that we introduce here in this latest of our series of Skyrme-HFB mass 
models, HFB-16, is the requirement that the contact pairing force reproduce exactly at each 
density encountered in the nuclear system (nucleus or the inner crust of neutron stars) the 
1 S'o pairing gap of neutron matter, as determined in microscopic calculations with realistic 
nucleon-nucleon forces. We retain the earlier constraints on the Skyrme force of reproducing 
the energy-density curve of neutron matter, and of having an isoscalar effective mass of 
0.8M in symmetric INM at the saturation density p ; we also keep the recently adopted 
device of dropping Coulomb exchange. Furthermore, the correction term for the spurious 
energy of collective motion has a form that is known to favour fission barriers that are in 
good agreement with experiment. 

Remarkably, despite the severe constraints imposed by neutron matter on both the 
Skyrme and pairing components of the effective force we have achieved the best fit ever 
to the mass data within the framework of mean-field models (the improvement is particu- 
larly striking for the most neutron-rich nuclei). The rms error a of the 2149 measured masses 
of nuclei with N and Z > 8 has been reduced to 0.632 MeV. Very few other mean-field mass 



tables have been published, the most successful of which is the RMF table of Ref. 16|, the 



rms error for essentially all measured nuclei being 2.1 MeV. Also to be noted is the mass 



table 
MeV 



aased on the SLy4 Skyrme force; this is limited to even-even nuclei and gives a = 5.1 
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39J. Turning to models lacking the microscopic basis of self-consistent mean- field 



models, we note that for the FRDM ("finite-range droplet model") 56| a = 0.656 MeV, 



while for the Duflo-Zuker model 



57l | we have a = 0.360 MeV (these last two values of a are 



both for 2149 nuclei). Only this last model does better than we do with the model presented 
here, but its applicability is limited to masses: it cannot be extended to any of the other 
quantities of astrophysical interest, such as fission barriers or the properties of the inner 
crust of neutron stars. 

Not only is our model well adapted to the determination of these quantities, it turns 
out that even with no flexibility at all remaining for the pairing force, the spectral pairing 
gaps that we find suggest that level densities in good agreement with experiment should be 
obtained. Finally, we find the correct sign for the isovector effective mass splitting without 
introducing additional terms in the effective force. 
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Much of the success of the present model must lie with the new constraint that we 
have imposed, i.e., our scheme for exactly matching the effective density-dependent contact 
pairing force to realistic microscopic calculations on neutron matter. We suspect also that a 
crucial role was played by the decision to make the pairing force for nucleons of charge type 
q depend only on p q . The modification of the collective correction (14 .2p has also contributed 
significantly to the improvement. (Dropping Coulomb exchange certainly accounts in part 
for the superiority of HFB-16 over HFB-14 and earlier models.) 

It is noteworthy that widely differing neutron-matter gaps can still correspond to the same 
finite-nucleus gaps, as well illustrated by the gaps for forces BSk7, BSk8 and BSkl6 shown in 
Fig.l, all of which lead to comparably good mass fits. This flexibility originates in the crucial 
role that the Skyrme force plays in passing from neutron-matter gaps to finite-nucleus gaps. 
It is this flexibility which allowed us to sucessfully impose the additional constraint of fitting 
to the neutron-matter gaps. On the other hand, it does not follow that one can start with 
any neutron- matter gap: for example, when we took screened neutron-matter gaps as our 
starting point we were unable to get good mass fits, no matter how the Skyrme force was 
adjusted. The fact that we do get better mass fits while neglecting medium effects suggests 
that these effects might be compensated in finite nuclei by, for example, coupling of particles 
via surface vibrations. Alternatively, the nuclear-matter medium effects themselves might 
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35]. 



be rather small, as indicated by recent quantum Monte Carlo calculations [32, 

To summarize, the HFB16 mass model not only gives a better fit to the mass data 
than any other mean-field model, but is also by far the most microscopically founded. It 
can thus be expected to make more reliable predictions of the highly neutron-rich nuclei 
of astrophysical interest. In particular, this is the first of our models well adapted to the 
investigation of a possible superfluid phase in the inner crust of neutron stars; in this respect 
it should be realized that since the model reproduces with precision all the mass data it 
must be giving a good account of surface properties, an important point in dealing with the 
inhomogeneities of the crust. 
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APPENDIX A: RELATION BETWEEN DIFFERENT FORMULATIONS OF 
THE HFB METHOD 

Computation of the HFB matrix elements. Our main concern in this appendix is to 
show how the computation of the oscillator matrix elements (12 . 1 1 ah and (|2.11bl) is greatly 
facilitated by starting with a density functional of the general form given by Eq. (12.161) . and 
working in coordinate space. First, however, we recall the definitions and various properties 
of the quantities appearing in this last equation. 

We have introduced 

(i) the nucleon density, 

<x=±l 

(ii) the kinetic-energy density (in units of h 2 /2M q ), 

T i( r )= 2 / d 3 r'5(r-r')V-V'p q (r,o-,r',a) (A.2) 

(iii) the spin-current density, 

J ff (r) = -i Y / d 3 r'S(r-r')Vp q (r,a;r',a')x<T a , a 

= i f d 3 r'6(r-T')V'p q {r,<T-y,o')x<T a , <T (A.3) 

a,a'=±l ^ 



and (iv) the so-called local pairing density [24j. |25] 



pi( r ) = p?( r ' a ; r > (7 )' ( A - 4 ) 



<T = ±l 



where <j aa > denotes the Pauli spin matrices, p(r,a;r' ,o~') and p q (r,o~;r,o~) are the normal 
and abnormal density matrices respectively, expressed in so-called coordinate space (position 
r and spin state a = ±1). From the general definitions of these two matrices 24, 25)], it 
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follows that the relationship between the coordinate space representation and the discrete- 
basis representation given in Eqs. (I2.7aj) and (I2.7bl) . respectively, is 

P q (r, a; r', a 1 ) = Pij ^(r, a)<f)*{r', a') (A.5a) 

and 

p q {r, a;r', a') = -a' re# <^(r, a)<f) j {r f , -a') , (A.5b) 

where 0j(r,cx) denote the s.p. basis wave functions. 

We return now to the question of the oscillator matrix elements ( I2.11al) and f!2.11bj) . 
Remarking that -Ehfb is a functional of p q (r) (and its gradient), r q (r), J q (r) and p q (r), 
which in turn depend on the matrices p^ and we have 

d -Ehfb 



dp 3 



d 3 r 



5 Erfb dp q {r) d£ RFB {r) dr q (r) <9£hfb(>) dJ q (r) 



5p q (r) dpji dr q (r) dpji dJ q (r) dp 



j 1 



(A.6a) 



and 



dE nFB _ f 3 <9£hfb0") <9p g (r) 



: _ 2££ = . 3 ^^v; ^ , (A.6b) 
JJ 54 J dp q (r) 8k* 3 1 J 

where we have introduced the functional derivative 

<5£hfb = d£ WB (r) _ d£ WB (r) 

5 Pq (r) ~ d Pq (r) ' dV Pq (r) ' 1 ' ' 

Substituting Eqs. (1 A. 1 1) - flA.3j) into Eq. (|A.6al) . together with Eqs. (IA.5al) . yields 



h V = E / "d 3 r^(r</)^(rU^-(ra) 

o-,cr'=±l ^ 



(A.8) 



in which h' q (r) aa / = aa'h' q (r)*_ al _ a is the self-consistent s.p Hamiltonian, appearing in the 
coordinate-space form of the HFB equations flA. 18|) 

h 2 

ti q {r) a , c = - V ■ 2M ^ (r) V 8 aa , + U q {r)6 aa , - iW q {r) ■ V x (A.9) 

where 

afe-W- "■m-Sct- ^"W- (A10) 
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To reduce Eq. (lA.6bl) in the same way requires that we make use of Eqs. ( ]A.4j) and 
( lA.5bl) . and also the expression 



(A.n; 



for the time-reversed conjugate &(r, a) of the s.p. state <fii(r, a). In this way we find 



A v = E / d 3 r0*(m)A ? (r)0 ? (m) 

<T=±1 J 



(A.12) 



in which A g (r) is the self-consistent local pairing potential, appearing in the coordinate-space 
form of the HFB equations (1A.18j) 

dS RFB (r) 



A ? (r) = 



dp q (r) 



(A.13) 



this must be real if time- reversibility holds, since p g (r) will then be real 2J, |25|. 

Explicit expressions for the fields appearing in Eqs. flA.lOj) and flA.13j) are as follows. 



h 2 dSnFB fi 2 



2M* 8r q 



1 

2M n + A* 1 



1 

+ - t2 



1 + \ Xl ) P ~ \\ + Xl J Pq 



(A. 14) 



and 



jj ^HFB ^ <9£hFB _ y <9£hFB 



9pg 



d(V Pq 



1 

+ J*' 



1 

+ -t 2 

4 



i \ / 3 . \ n 



1 + ^i ) ( T ~ o V P ) - I 2 4 ' 1 ) I r ' 



2 
1, 



" 2 
1 



+ 12 t3 



(l + ^s) ( 2 + 7)P 7+1 - Q + *s) (vp« + 7P 7 ' 1 E pj 

Vo(v-j + v.j 9 ) + ^ c - + i ^ ^ 



4 ^ 

q—n,p * 



W q = = -W V(p + p ? ) - -(txX! + t 2 X 2 )J + -(ti - t 2 )J, 



(A.15) 



(A.16) 
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Note that for the Coulomb field term S qtP V Coul appearing in Eq. (1A.15I) we have taken the 
electrostatic potential given by Eq. (12.181) . which, in view of the finite proton size, is not 
given exactly by the functional derivative of the Coulomb energy (5Ec OVL \)/(Sp p ). 
The local pairing field, defined in Eq. flA.131) . is given by 

\ = ^ = ^[ Pn ,p p }p q . (A.17) 

Let us point out that the pairing functional £ pa i r contributes not only to this pairing field 
but also, through its dependence on the nucleon density, to the s.p. Hamiltonian h' (r) aa > 
as the last term of Eq. (1A.15j) . which is essentially a rearrangement term, corresponding to 
the second term of Eq. (I2.14p . (All the other rearrangement field terms, corresponding to 
the first term of Eq. (I2.14p . appear in Eq. (1A.15I) . and can be identified by their coefficient 

hi.) 

Transformation from discrete-basis to coordinate-space formulations of HFB equations. 
It is convenient now to point out that the HFB equations in coordinate space can be ob- 



tained directly from the standard formulation of Refs. [19|, [20J], i.e., from the corresponding 
equations in discrete-basis form (I2.10p . by substituting in the latter Eqs. (1A.8I) and (1A.12j) . 
and invoking time-reversibility. 

This leads for each nucleon species q to the coordinate-space form of these equations, as 



given by Refs. 



24 



25|], thus 



A q (r)5 aal -h> q {r) a „, + \ q 8^)\^{r,a>)) \1${r,&)) ' 

where the s.p. Hamiltonian h q (r) aa > and pairing field A 9 (r) for the forces (12. ip and (12.21) are 
as given in Appendix |A] (note that the rearrangement terms are included here). The upper 
and lower components of the Bogoliubov quasiparticle wave function, denoted respectively 
by ip[f(r, a) and ip^i{r, a), are related to the U and V matrices introduced in Eq. (12. 5p by 

4t (r,a) = J2 Up ^ (r, a) , 4f (r , a) = ^ V Jt ^ (r , a) ; (A.19) 

Hi) 

note the appearance of the time-reversed conjugate state in the latter expression (Eq. (6) of 
Ref 58|] is wrong). 

Using Eqs. ( 1A.19I) . (1A.5ap . (1A.5bl) and the assumed time-reversal invariance, we have the 
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identities 



Pg 



(A.20) 



»(«) 



and 



p q (r,a;r',a') = - £ ^ (r, a)^f(r', a')* = - J>i?(r, a)V&V, <r^ 



(A.21) 



Applying the properties of the U and matrices, as given for example by Eqs. (7.5) of 



Ref. 



201 ]. to the definition (IA.19I) it can be easily checked that the quasiparticle wavefunction 



satisfies the completeness relations given by Eqs. (2.20a) and (2.20b) of Ref. [24]. 



APPENDIX B: HFB EQUATIONS IN UNIFORM MATTER 

In a uniform system, it is natural to replace the oscillator basis by a plane-wave basis, 

(fiki?, <j) = \= exp (ik ■ r) x( cr ) > (B.l) 
V v 

where %{p) is the Pauli spinor and V is the normalization volume. It is easily seen from 
Eq. (1A.8|) that the s.p. Hamiltonian h' (r) aa i is diagonal in this basis, 

h' kl = 4 ?) Ski , (B.2) 

with 

Also, since the pairing field A 9 (r) must be independent of r in a uniform system it follows 
from Eq. (1A.12I) that all matrix elements A/& vanish unless I = k, i.e., 

A lk = 5 fk A q = -A kl (B.4) 

where A q is a constant for a given uniform system at a given density, given by Eq. (1A.17|) . 
The assumption that the pairing force acts only between states of the s.p. Hamiltonian that 
are the time-reversal of each other constitutes the essence of the BCS method (note that 



this approximation differs from that discussed in Ref. 25|), an approximation to the HFB 



method that is often convenient for finite nuclei, but we have shown here that for uniform 
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systems the HFB method reduces exactly to the BCS method. Furthermore, for infinite 
but inhomogeneous nuclear matter, the weaker the departure from homogeneity the better 
the BCS method will approximate the HFB method; this means that in the inner crust of 
neutron stars the BCS approximation will be better the deeper the layer in question. 

More precisely, by inspecting Eq. f]A.12|) and by remembering that only s.p. states close 
to the Fermi level contribute to pairing correlations, it can be seen that in general the BCS 
approach is justified whenever the pairing field is slowly varying in the spatial domain for 
which the s.p. wavefunctions for states around the Fermi level, are non-vanishing. This 
condition is usually fulfilled for strongly bound nuclei (for which the chemical potential X q 
lies deep inside the s.p. potential well) since the pairing field A g (r) is typically more or less 
constant in the nuclear interior, where the bound s.p. wavefunctions take their largest values. 
In contrast, for weakly bound nuclei, pairing correlations involve not only bound states but 
also states from the continuum. Since the pairing field vanishes outside the nucleus, it follows 
that it will be varying rapidly in a region where the s.p. states that contribute to pairing are 
still strong, whence the BCS approximation can be expected to break down for such nuclei. 

Solving now the HFB equations ( I2.10j) by using Eqs. ( IB. 21) and ( IB. 41) yields 



4" = v/(4'' - \) 2 + ^ , (B.5a) 



Ult = U$ = -L + fl^A ] (B.5b) 

and 




1/2 



making use of the properties of the U and V matrices. Using Eq. ( 1A.19I) . ( 1B.5bl) and ( lB.5cl) 



it can be seen that the quasiparticle wavefunction in uniform nuclear matter reduces to 

^ } (r, a) = Mr, <r) , 4k(r, *) = V u Mr, a) , (B.6) 

where (pk( r , o") is given by Eq. (IB. 1|) . 

Substituting Eqs. (lB.5b|) and (lB.5c|) in Eqs. (12.7aj) and (I2.7bl) leads to the familiar ex- 
pressions of the normal and abnormal density matrices 

Pu = (VkkY ^ki , Kki = V k iU k k hi ■ (B.7) 
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Using Eq. (1A. 1[) and (1A.5aj) together with (IB. II) and (1B.5cl) . the nucleon density p q (r) = p q 
is given by 



+ /_ £ _ Xq y 



Pi" 4vr 2 



Likewise from Eqs. ( 1B.7I) and ( I2.9bl) we find the "BCS gap equation" 



1 f 2M ;\ 3/ \.n 0t „ „ 1A /\ 



n 2 \ h 



A * = KtM ^1Pn,P P ] A, / d £ - v (B.9) 



> " \? + A^ 



The integral is taken inside the subspace A introduced to regularize the ultra-violet diver- 
gences, arising from the zero range of the pairing interaction. These divergences are removed 
by imposing a cutoff, either in the q.p. energy spectrum or in the s.p. energy spectrum. In 
the latter case, different prescriptions have been employed in the literature: (i) e < E\, (ii) 
£ < U q + E\, (iii) £ < X q + E\, (iv) \ q — e A < e < \ q + where £\ is a constant. The fixed 
cutoff (i) is the only choice which implies that the gap equations depend explicitly on the 
s.p. potential U q . Since a density dependent pairing force affects directly the s.p. energies 
via the rearrangement potential (see Eq. (1A.15|) ). the gap equations (1B.9|) will involve not 
only the pairing strength v nq [p n ,p p ] but also its partial derivatives. Note that since the 
integrand in the integral appearing in Eq. (1B.9j) takes significant values only in the vicinity 
of the Fermi level, the cutoff prescriptions (i) and (ii) lead to essentially zero pairing gaps 
A q ~ for densities such that ea < X q and £\ < p q respectively. A peculiar consequence 
is that the pairing gap may be vanishingly small while the pairing force is not. We have 



encountered such a situation for the Skyrme force SLy4 with the pairing force of Ref. 
using the cutoff prescription (i), as can be seen by comparing Figs. [[]and[2j 
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TABLE I: Parameters of the analytical fit (|3.2p for the BCS pairing gap in neutron matter extracted 
from Fig. 7 of Ref. 28j (the unit of length is fermi and the unit of energy is MeV). 



A 


h 


k 2 


^3 


h 

ft max 


910.603 


1.38297 


1.57068 


0.905237 


1.57 



TABLE II: Force BSkl6: lines 1-10 show the Skyrme parameters, lines 11-14 the pairing parameters 
and the last four lines the Wigner parameters (see Section HVl for further details). 



uq iivie v nn j 


-1 837 93 

IOO 1 . 


1 1 riVTpV fm 5 l 

u J_ 1V1G V 1111 


383 521 


62 uviev im j 


q 41 7"\P> 


63 [iviev tm 'J 


1 1 ^9^ n 


Xo 


n a "^ofinn 

U.4loZDUU 


X\ 


D 8941 OR 


X2 


44.D0ZU 


X3 


0.689797 


Wq [MeV fm 5 ] 


141.100 


7 


0.3 


fn 


1.06 


ft 


0.99 


fp 


1.05 


E\ [MeV] 


16.0 


Vw [MeV] 


-2.60 


A 


240 


V w [ MeV l 


0.70 


A 
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TABLE III: Parameters of collective correction for model HFB-16 (see Section IIVI for further 
details). 



b (MeV) 


0.8 


c 


10 


d (MeV) 


2.6 


I 


10 


ft 


0.1 



TABLE IV: Rms (a) and mean (e) deviations between data and predictions for model HFB-16; 

n n 

for convenience we also show models HFB-15 [111 ]. HFB-14 [10( and HFB-8 [7|. The first pair of 
lines refers to all the 2149 measured masses M, the second pair to the masses M nr of the subset of 
185 neutron-rich nuclei with S n < 5.0 MeV, the third pair to the neutron separation energies S n 
(1988 measured values), the fourth pair to beta-decay energies Qp (1868 measured values) and the 
fifth pair to charge radii (782 measured values). The last line shows the calculated neutron-skin 
thickness of 208 Pb for these models. 





HFB-16 HFB-15 HFB-14 HFB-8 


a(M) [MeV] 


0.632 


0.678 


0.729 


0.635 


e(M) [MeV] 


-0.001 


0.026 


-0.057 


0.009 


a(M nr ) [MeV] 


0.748 


0.809 


0.833 


0.838 


e(M nr ) [MeV] 


0.161 


0.173 


0.261 


-0.025 


a(S n ) [MeV] 


0.500 


0.588 


0.640 


0.564 


e(S n ) [MeV] 


-0.012 


-0.004 


-0.002 


0.013 


a(Qp) [MeV] 


0.559 


0.693 


0.754 


0.704 


KQp) [ MeV i 


0.031 


0.024 


0.008 


-0.027 


cr(R c ) [fm] 


0.0313 


0.0302 


0.0309 


0.0275 


e(R c ) [fm] 


-0.0149 


-0.0108 


-0.0117 


0.0025 


^( 208 Pb) [fm] 


0.15 


0.15 


0.16 


0.12 
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TABLE V: Macroscopic parameters for force BSkl6 (for convenience we also show forces BSk 
15 11] and BSkl4 [10j] ) . The first twelve lines refer to infinite nuclear matter, the last two to 
semi-infinite nuclear matter. See Section [IV] for further details. 





BSkl6 


BSkl5 


BSkl4 


a v [MeV] 


-16.053 


-16.037 


-15.853 


PO [fm~ 3 ] 


0.1586 


0.1589 


0.1586 


J [MeV] 


30.0 


30.0 


30.0 


M*/M 


0.80 


0.80 


0.80 


M*/M 


0.78 


0.77 


0.78 


K v [MeV] 


241.6 


241.5 


239.3 


L [MeV] 


34.87 


33.60 


43.91 


Go 


-0.65 


-0.67 


-0.63 


G 


0.51 


0.54 


0.51 


Gi 


1.52 


1.47 


1.49 


G\ 


0.44 


0.41 


0.44 


Pfrmg/PO 


1.24 


1.24 


1.24 


a s f [MeV] 


17.8 


17.7 


17.6 


Q [MeV] 


39.0 


39.7 


35.0 



TABLE VI: Rms and mean differences between predictions for highly neutron-rich nuclei (4.0 MeV 
> S n > and 26 < Z < 110) given by different pairs of mass models. Mean differences are shown 
in parentheses. 





M S n Q p 


HFB-15 - HFB-16 
HFB-14 - HFB-16 
HFB-8 - HFB-16 


1.356 (0.977) 0.378 (-0.033) 0.520 (0.092) 
3.230(-2.465) 0.392 (0.136) 0.595 (-0.343) 
1.581 (0.319) 0.581 (-0.099) 0.892 (0.198) 
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